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We present a variational study of pseudo-spin 1/2 Bose gases in a harmonic trap with weak 
3D spin-orbit coupling of <x ■ p type. This spin-orbit coupling mixes states with different parities, 
which inspires us to approximate the single particle state with the eigenstates of the total angular 
momentum, i.e. superposition of harmonic s-wave and p- wave states. As the time reversal symmetry 
is protected by two-body interaction, we set the variational order parameter as the combination of 
two mutually time reversal symmetric eigenstates of the total angular momentum. The variational 
results essentially reproduce the 3D skyrmion-like ground state recently identified by Kawakami et 
al. We show that these skyrmion-like ground states emerging in this model are primarily caused by 
p wave spatial mode involving in the variational order parameter that drives two spin components 
spatially separated. We find the ground state of this system falls into two phases with different 
density distribution symmetries depending on the relative magnitude of intraspecies and interspecies 
interaction: Phase I has parity symmetric and axisymmetric density distributions, while Phase If 
is featured with special joint symmetries of discrete rotational and time reversal symmetry. With 
the increasing interaction strength the transition occurs between two phases with distinct density 
distributions, while the topological 3D skyrmion-like spin texture is symmetry protected. 

PACS numbers: 67.85.Fg, 03.75.Mn, 67.85.Jk, 03.75.Lm 


I. INTRODUCTION 

The experimental realization 0,1 of one-dimensional 
(ID) spin-orbit (SO) coupling in pseudo-spin 1/2 Bose 
gases has stimulated many theoretical works on SO cou¬ 
pling in cold atom physics. These works range from Ra¬ 
man induced ID SO coupling that has been real¬ 

ized in cold atoms to more symmetric two-dimensional 
(2D) Rashba configuration |7Hl5j that has been exten¬ 
sively studied in condensed matter. In the absence of 
harmonic trap, single particle ground states of both Ra¬ 
man induced and Rashba SO coupling are degenerate, 
and two-body interaction selects the generic ground state 
from the degenerate manifold determined by the interac¬ 
tion parameters. For example, Wang et al. Q found 
two distinct ground state phases, namely the plane wave 
and standing wave(or stripe) phases, appeared when in¬ 
traspecies two-body interaction is larger or smaller than 
interspecies interaction respectively in homogeneous 2D 
Rashba SO coupled pseudo-spin 1/2 Bose gases. In the 
presence of a 2D harmonic trap, a more complex phase di¬ 
agram of Rashba SO coupled Bose gases with two classes 
of phases and several subphases in each was figured out 
by Hu et al. 00- 

Now experimental schemes for the realization of 
Rashba SO coupling have been proposed such as in Ref. 
0. On the other hand, the most symmetric three- 
dimensional (3D) SO coupling or Weyl coupling, which 
even doesn’t exist in solid matter, is expected to be re¬ 
alizable in cold atoms gases and experimental schemes 
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for that have also been proposed theoretically 0i- 
Recently, Kawakami et al. [19j identified a 3D skyrmion 
ground state in 3D SO coupled two-component bosons 
by numerically minimizing the Gross-Pitaevskii (GP) en¬ 
ergy functional of the system. They explained the sta¬ 
bility of 3D skyrmion ground state as a result of helical 
modulation of the order parameter in the presence of 
SO coupling. The interaction in their work is supposed 
to be SU(2) symmetric. Even though two skyrmion-like 
ground states are found to be stabilized in different in¬ 
teraction regimes, a ground state phase diagram is still 
absent now. In another work by Li et al. 20], the 3D 
skyrmion-like ground state is found to emerge in weak SO 
coupling regime, while skyrmion lattice arises in strong 
SO coupling regime. 

In this work, we consider a pseudo-spin 1/2 boson sys¬ 
tem subject to 3D SO coupling of er ■ p type in a har¬ 
monic trap, and aim to elucidate the role of interaction 
in determining the ground state density and spin texture 
therein. In Sec. II we introduce the energy functional for 
the model in rescaled units of length, energy, interaction 
and SO coupling strength. In weak SO coupling case, 
the single particle energy levels are essentially harmonic 
oscillator-like 0l|, and SO coupling will mix states with 
different parities while keeping the total angular momen¬ 
tum a conservative. In Sec III we first try to couple two 
lowest s and p wave states with the same total angular 
momentum 1/2 into two spinor wave functions with to¬ 
tal angular momentum magnetic quantum number ±1/2 
that are time reversal state of each other. Then we set 
the variational order parameter as superposition of these 
two states just as has been done in ID and 2D cases 
[3, 0]. Finally, we calculate the energy functional us¬ 
ing the proposed variational order parameter. In Sec IV 
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the ground state phase diagram is determined by numer¬ 
ically minimizing the energy functional with respect to 
the variational parameters, we illustrate the density and 
spin texture for the two phases. Sec. V summarizes our 
main results. 


II. MODEL 


We consider a pseudo-spin 1/2 boson system confined 
in a harmonic trap with a weak Weyl type 3D spin-orbit 
(SO) coupling er-p. The system is described by its Gross- 
Pitaevskii (GP) energy functional under the mean-field 
approximation 


£ — £q + £int , 

where the single particle part is 


( 1 ) 


£o = J d 3 (r) — f -muj 2 r 2 + Xcr ■ T (r) (2) 


with m the mass of atoms and ui the trap frequency. 
T = (ip-\,ip\) T denotes spinor order parameters for 
bosons with pseudo spin states <r = (a x ,a y ,a z ) are 
the Pauli matrices and A parameterizes the SO coupling 
strength. The interaction takes the usual contact 
form of s-wave scattering interaction j23j. We assume 
now [M, [13 the two intraspecies interaction parame¬ 
ters being the same g^ = = g and define the relative 

magnitude of the interspecies and intraspecies parame¬ 
ters as c = g'l'l/fftT- The interaction part is then 

£int = ^ J d 3 r ((g + eg) n 2 +4 (g - eg) S 2 ) . (3) 


In Eq. ([3]), n (r) = n^ (r) + n± (r) is the particle den¬ 
sity and S z is the 2 component of the spin density 
S = i'pl'er'I' with (r) = It/ti ( r )| 2 the particle den¬ 
sities of two components, respectively. The correspond¬ 
ing Hamiltonian is time-reversal (TR) symmetric with 
time reversal operator defined as T = — icr y K and K 
denotes the complex conjugate. The system has length 
scale of the trapping potential It = \/h/moj , energy scale 
hco, interaction strength scale Hcol^/N, and SO coupling 
strength scale \Jhui/m. If we further normalize the order 
parameter to unity, i.e., 4^ —> with N the total 

particle number in the condensate, the energy functional 
per particle is obtained as 

e = J d 3 r& (r)|-^- + ^- + A(T -p|^( r ) 

+ i / d3r ^ 9 + C ^ ^ + 4 ^ ~~ ° 9>> S *) ' ^ 


III. VARIATIONAL APPROACH 

In the case of weak SO coupling the single parti¬ 
cle energy spectrum in our system should be harmonic 


oscillator-like as proposed in Ref. J2(| |2l|. The three 
dimensional harmonic oscillator thus proves to be a good 
choice of the trial wave function, upon which we may 
develop our variational method. As can be seen later 
the spin-orbit coupling induces transition between eigen 
states with the same total angular momentum but dif¬ 
ferent parity, which are mixed into the variational wave 
function. The interaction Hamiltonian further couples 
the two time-reversal states with different weight factor 
due to the anisotropic interaction parameter ratio c. 


A. Variational order parameter 

The eigenequation of three dimensional harmonic os¬ 
cillator, Y —K 4> = £</>, has well-known solutions, 

with energy eigenvalues e„ r ; = 2 n r + l + | and eigen¬ 
functions (pn r im, (r,0,tp) = Rn r i (r)Yi mi (6, ip). Here n r 
is the radial quantum number, l is the orbital angular 
momentum quantum number with mi its magnetic quan¬ 
tum number, R Ur i is the radial wave function, and Yt rni 
is the spherical harmonics. The Casimir operator l 2 and 
s 2 for the orbital and spin angular momenta and their 
2 -components are all conservatives in the harmonic os¬ 
cillator problem. In order to take into account the spin- 
orbital coupling term <x • p, it is convenient to choose the 
coupled representation of angular momentum, i.e. the 
complete set of commutative operators l 2 , s 2 , j 2 , j z where 
j = 1 + s and j z denote the total angular momentum 
and its 2 -component, respectively. The eigenfunction 
< j>n T imi ( r , 0, f) should be combined with the spin wave 
function Xm„ in the coupled representation as 

(pn r ij mj ( r , 0, tp) = R nr i (r) Yj m . (H), (5) 

where Yj (fi) = Y lmi Xm s is the spinor 

spherical harmonics [24] with j = l± 1/2 and C^ m \ the 

Imi -^rris 

Clebsch-Gordan coefficients. In the coupled representa¬ 
tion, the ground state wave function has n r = l = 0. This 
gives a total angular momentum j = \ with rrij = 
and the two degenerate ground states are 

</>oo±±± ( r ) = R oo ( r ) y| ± i (H), (6) 

respectively. Because the SO coupling term breaks the 
parity symmetry, it can couple s and p wave states with 
the same total angular momentum j and j z [20j . Keep¬ 
ing these consideration in mind, in the simplest approx¬ 
imation, we suppose the ground state contains only the 
lowest s and p wave states with total angular momentum 
quantum number j = A in presence of the SO coupling 
term. The state with mj = A takes the form 

^=i, raj = i = N <* (</>00±± + W>01±±) ( 7 ) 

where N a = (1 + a 2 ) -1 / 2 , a stands for the relative weight 
of the s and p orbital modes, and i in front of a originates 
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from the pure imaginary matrix element of the SO cou¬ 
pling between the two states in Eq. 0. This hypothesis 
is similar to that appears in Refs. (?() and , and has 
been verified numerically (2(j. Explicitly this state is a 
spinor 


$ = i i = N a 

J 2 ' 3 2 


RqqYqo — ioi\/\R q\Yw 


iot\l o-Rqi^ii J 


( 8 ) 


The state with rrij = — i takes the form 


$j == N <* V‘/'oo|-A 


i _i + iad> m i _ i 

o o ' JJ - o o 


) (9) 


and similarly we have 


_i m __i = N a 


—iay 3-RoiTi-i 
Rookbo + 4i?oikio 


( 10 ) 


which is nothing but the time reversal of -_i _i. In 

J — 2 ’ " l 3 ~ 2 

the single particle level, i jTOj= ±i and any normal¬ 
ized superposition of them has the same energy thus are 
“degenerate” single particle states, which is similar to 
degeneracy indicated by Kramers’ theorem in spin-1/2 
system. 

The single particle states exhibit infinite-fold degen¬ 
eracy and we expect this degeneracy can be partially 
resolved by the interaction which would pick up the 
ground state from these degenerate states as in the case 
of Rashba spin-orbital coupling considered by Wang and 
Zhai Q. Since the interaction doesn’t break the time 
reversal symmetry, the residual two-fold Kramers degen¬ 
eracy need to be considered in the wave function [1]. We 
therefore set the variational order parameter as 


T = c + 4> + c_T$ 

= f c + $ t-c-$n 
V C + ^ + c_$| J ’ 


( 11 ) 


with the constraint cl + c 2 =1. Here $ = <K_ i m _i 

i J — 2 ,r,L 3 — ~2 

and are its up and down components. So far, we 
have introduced three variational parameters a,c+,c_ 
and the energy functional of Eq. 0 can be calculated 
analytically using the proposed order parameter CO- 


B. Energy functional 

We calculate the energy functional on the variational 
wave function CO- The contribution comes from two 
parts, the single particle and the interaction Hamilto¬ 
nian. We notice that for the kinetic and trapping poten¬ 
tial terms the nonzero integral contribution comes from 
those states with the same parities, while spin orbital 


coupling er ■ p term will mix states with opposite pari¬ 
ties, i.e. 



Here we have used 


^OQii ' Pi '/’Olii 



which is on account of \j z , cr ■ p] =0. 

It is crucial to calculate the contribution of the spin- 
orbital coup ling term by means of the irreducible ten¬ 
sor method [24j . To this end we first introduce the ir¬ 
reducible form of spin-orbital coupling term. The irre¬ 
ducible tensor form of momentum operator is (24 ] 

p^ = iV2^{c^l^} (1) (14) 

where C ^ and /W are rank-1 irreducible tensors of the 
unit vector r and the orbital angular momentum 1, and 

defines the rank-fc tensor product of rank- 
m irreducible tensor A and rank -n irreducible tensor 
B^ n \ According to jUj], the dot product of two arbitrary 
vectors A and B is related to the tensor product through 
A ■ B = — \/3 {A.W.BW}*' 0 ). In our case, the radial co¬ 
ordinate r can be separated from the spin and spherical 
parts accordingly 


tr p = /„<!) jcWi* 1 *} 11 ’!. 


(in 


( 0 ) 


(15) 


such that 


^ooii I* 7- ‘ Pi roiii 


iV& (^Rqq (r) 


x ( Y? k (Cl) 
' 2 2 


r« 


Roi ( r ) 
jcKD/Mj 


"} 


(0) 


Fi 1 ! (Cl) 


+ iV 3 (^Roo (r) 
x (y ?i (Cl) |{cr (1) C ,(1) } 


Roi ( r ) 
( 0 ) 


Yh (Cl)). 
22 ' 


( 16 ) 
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The integrals for the radial coordinates are easy to cal¬ 
culate, 


Roo (r) 


Roo ( r) 


1 

r 

d_ 

dr 


Roi ( r ) ) = 


R 11 (r)) = - 7S , 


(17) 


(18) 


where Roo (r) = \/2 2 /y/ne r / 2 and R 0 1 (7') = 
•\/2 3 / (3y / 7r)re _r / 2 are used. Wigner-Eckart theorem 
can be used to calculate the angular and spin integral 


Y ?i {n) 
2 2 


(DjcW/W} 


(i) 


( 0 ) 




(19) 



FIG. 1: Phase diagram of weakly SO coupled two-component 
Bosons with coupling strength A = 0.2, which shows two 
skyrmion-like phases I and II. Phase I is a skyrmion ground 
state of order parameter exp [— id (r) • S] £ with £ 2 = (1, 0) T 
and Phase II is a skyrmion state with £ x = -^=(1, 1) T - Density 
distribution and spin texture of these two phases are shown 
in Fig. [2] and Fig. [3] respectively. 


Y? i (ft) 




(o) 


Yh 


(!! > > = 7r 


( 20 ) 


Substitute Eq. (11711201) into Eq. (TTB1) . one has 


^oo ii 




( 21 ) 


Hence the single particle part of the energy functional is 
J dW (r) + T — + Act • p J 'P (r) 

= K 2 + VeaX^j ( 22 ) 


Collecting Eqs. (l22l) . (l24l) and (l26l) into Eq. Q, we 
finally arrive at the variational result for the ground state 
energy per particle 

e = + -a 2 + v^aA^j 

3 

+ N*( 2 n)- 2 — [(11 a 4 + 12a 2 + 36 )g + (4 a 4 + 24 a 2 )cg 
(g — eg) ( a 4 — 12cx 2 + 12)] . 

(27) 


where we have used the eigen energies of the s and p 
states of the three dimensional oscillator are respectively 
£oo = 3/2 and eoi = 5/2. 

For the calculation of the interaction part of energy 
functional, it is easy to show that the total density n is 
always spherical symmetric 

n = |$| 2 = (d^)- 1 ^ 2 (i? 2 0 + a 2 Roi) , (23) 

and the density-density interaction energy is 

J d 3 rn 2 = 27r)"i(5a 4 + 12a 2 + 12). (24) 

On the other hand, the spin density is anisotropic, e.g. 
the 2 component takes the form of 

S z = (87r) _1 iV 2 {(c+ - c 2 )(i?Q 0 + a 2 R 2 01 cos2 6) 

— (2c + C-)(2aRooRoi sindsin^ — a 2 R$ 1 sin 29 cos ip) } , 

(25) 

and the spin-spin interaction energy is integrated as 

J d 3 rS 2 = iA^(27r)“i [^(7a 4 - 12a 2 + 36) 

-\c 2 + c 2 _{a 4 - 12a 2 + 12)]. (26) 


IV. GROUND STATE PHASE DIAGRAM 

The ground state phase diagram can be determined 
numerically via the minimization of the variational en¬ 
ergy with respect to the parameters a, c+ and c_ for 
given c and g. We notice that the parameters c+ and c_ 
appear only in the last term of Eq. m in the form of 
c\c 2 _, the value of which ranges from 0 to 1/4. The 
parameter c^_c?_ as a whole takes the value of either 
0 or 1/4 in the minimization, depending on the signs 
of ( g — eg) and /(a) = a 4 — 12a 2 + 12. The ground 
state thus falls into two classes of phases as depicted in 
FIG-IU Phase I, the variation yields |c+| 2 = 1, |c_| 2 = 0 
or |c + | 2 = 0, |c_| 2 = 1; Phase II, the variation yields 
|c+| 2 = |c_| 2 = 1/2. It is clear that a must be negative 
for a positive A due to the fact that a’s in Eq. (EH) are 
all even-ordered except the spin-orbital coupling term. 
We see that c = 1 divides the phase plane into upper 
and lower regions. With increasing g the system enters 
alternately into Phases I and II and the boundaries are 
determined by /(a) = 0, i.e. a± = — \/6 ± 2\/6. For 
typical experiments with bl Rb condensate, the interac¬ 
tion strength scale is 10 _13 Hzcm 3 , which gives rise to 
g ~ 40 — 80. For smaller g, a > a_ results in a pos¬ 
itive value of /(a) such that c > 1 region belongs to 
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time reversal states of each other and have similar density 
and spin texture except that the spin-up and spin-down 
components are exchanged. The order parameter for the 
former has the form 


= (47t) 2 N a 


Roo (r) — iaRoi (r) cos 9 
—iaRoi (r) sin 9e l ' p 


The particle densities for the spin-up and spin-down com¬ 
ponents are 

n t = ( 47 r ) —1 (Rio (r) + a 2 R 2 01 (r) cos 2 6) , 
n i = N 2 a 2 R^ (r) sin 2 0, (28) 


which respect the rotational symmetry about z-axis and 
the parity symmetry, i.e. 71 ^, 4 . (r) = n^^(r,6) and 
6) = n^ t i(r, tt — 9). The densities of the two com¬ 
ponents in xz and yz planes are the same as shown in 
Fig. 1 which exhibit clearly characters of the p wave 
state, i.e. the spin-up component is dumbbell-like while 
the spin-down component forms a torus. The total den¬ 
sity on the other hand is isotropic - the sum of n-j- and 
714 in Eq. (1281) relies only on the radius r. 

This spin density calculated on the variational order 
parameter shows interesting spin texture described by 

S x = ( 47 t) ~ 1 N 2 (aR qo (r) R 01 (r) sin 9 sin ip 
+a 2 i?oi ( r ) s i n ^ cos 0 cos p) j 
S y = ( 47 r) _ 1 iV 2 (—ai?oo (r) R 01 (r) sin 9 cos ip 
+a 2 i?o! (r) sin 9 cos 9 sin ip) , 

S z = ( 87 r) _ 1 iV 2 (i?Q 0 (r) + a 2 !?^ (r) cos 29). (29) 


FIG. 2: (Color online). Density distribution and spin texture 
of Phase I for a = —1.6. Top: Three rows are densities in 
xy, yz and xz planes respectively, three columns are for up, 
down components and the total density as explicitly labeled 
above each column. Density distributions in xz and yz plane 
are the same due to the z-axis rotational symmetry. Bottom: 
3D skyrmion spin texture s(r) = S(r)/n(r) of Phase I. The 
streamline plot of s in a selected region are shown. 

Phase I and c < 1 belongs to Phase II. By adjusting the 
trapping frequency and the density of the condensate one 
can easily increase g to cross the critical line such that 
a £ [a+,a_], which makes f[a) negative. We observe 
an interesting swap of the phases: c < 1 corresponds 
to Phase I and c > 1 corresponds to Phase II. Similar 
phase transition appears in the Rashba spin-orbital cou¬ 
pled Bosons [9|, 05 • Further increasing the interaction 
strength makes the optimized parameter a < a + and 
the phases swap occurs again. 

The density distributions and the spin texture of 
Phases I and II are shown in Figs. [2] and [3] respectively. 
Typical features include: 

Phase I: this phase contains two degenerate states 
|c+| 2 = 1, |c_| 2 = 0 and |c + | 2 = 0, |c_| 2 = 1. They are 


The average value of the spin in the xy plane is zero, i.e. 
{S x ) = (S y ) = 0. The spin texture s(r) = S(r)/n(r) is 
depicted in Fig. [2] and we find that spin density forms 
a torus near the xy plane and a bundle of nearly ver¬ 
tical streamlines of spin penetrate the central region of 
the torus. This skyrmion-like texture has been discussed 
in Ref. 0 and identified as the ground state in c < 1 
regime for an interaction parameter cq = 100. Li et al. 
[20( | also found this ground state skyrmion spin texture 
in weak SO coupling case for isotropic interaction c = 1. 
The term “skyrmion-like” means the absence of bound¬ 
ary condition at r —> 00 [19] thus the winding number 
for the texture is not an integer. 

In order to get a deep understanding of the skyrmion 
nature of this ground state, we notice that the order 
parameter can be obtained from a local spin rotation 
from the polarized spinor wavefunction = (c+, c_) T = 
(1,0) T 

T- = exp {-in (r) -s) sjn (r)£z (30) 

supposing that n (r) = (dTr^IV 2 (l?g 0 (r) + a 2 I?oi ( r )) 
and ft (r) = u>(r)r/r. This operation rotates the spin 
at position r by an angle u (r) about the axis r/r. 
The rotation angle is position dependent, i.e. ui (r) = 
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2 arctan (ai?oi i r ) /Roo ( r )) and s is the usual spin an¬ 
gular momentum operators for spin-1/2. It is the ex¬ 
plicit form of ft (r) that determines the specific texture 
of the skyrmion [221 ]. The polarized spinor order pa¬ 
rameter C, z has all spins being oriented in the positive 
z-direction. After the rotation the order parameter Tc 
for this skyrmion state is position-dependent. The order 
parameter is the most symmetrically shaped skyrmion 
with the symmetric axis unrotated, which is identical 
with that already discussed in Refs. 0 ii-0. In 
those papers the 3D skyrmion states are proposed as ex¬ 
cited states in pseudo-spin 1/2 or ferromagnetic spin-1 
0 Bose gases although they may be metastable. 

Phase II: this phase again contains two degenerate 
states with c+ = c_ = ±-^ and c. (_ = — c_ = ±^=, 
which are time reversal states of each other. They share 
similar density distribution and spin texture just like the 
case of Phase I. The order parameter for c+ = c_ has the 
form 


T = (87r)-^TV a 


R 00 (r) - iaRoi (r) (cos 9 + sin de 1<P ) \ 
Roo ( r) — iaR 0 i (r) (— cos 9 + si\\9e l ' p ) J 

(31) 


and the particle densities for the two components are 

n t = ( 8 tt)~ 1 N 2 (R 2 0 (r) + a 2 R 2 01 (r) 

—2aRooRoi sin 9 sin p + a 2 Rh sin 2 9 cos p) , 

H = ( 87r )“ l7V « (#00 ( r) + a 2 R^ (r) 

+2aRooRoi sin 9 sin p — a 2 Rh sin 29 cos p) . (32) 

The density for each component consists of two parts, 
one is isotropic that is common for both components, 
the other is complementary to each other as shown in 
the second and fourth lines in Eq. (1321) . This leads 
again to an isotropic total density. The overall density 
distribution of the two components can be visualized as 
two cashew nuts perpendicularly crossing and partially 
overlapping with each other. The distributions in xy , yz 
and xz planes are shown in Fig. [3] The density dis¬ 
tributions have the following symmetries, e.g., the den¬ 
sities of two components are invariant under the com¬ 
bined operation of time reversal and tt rotation about 
a;(or z) axis, i.e. n^(r, n — 9,2n — <fi) = ni(r,9,<f>), 
Tif(r, 9 1 7T + cp) = nj,(r, 9 , </), while the n rotation about y 
axis itself leaves the density distributions unchanged, i.e. 
%,!(D tt - 9, 7T - 4 >) = n t)4 .(r, 9 , </). 

The spin texture associated with the order parameter 
is expressed as 

S, = (Stt)" 1 !^ 

(l?g 0 + a 2 Roi (sin 2 0cos2y> — cos 2 9 )) , 

S y = (8n)- 1 N 2 x 

(2aRooRm cos 9 + sin 2 9 sin 2p) , 

S z = (87t) -1 1V 2 x 

(— 2aRooRo\ sin 9 sin p + a 2 R^ ± sin 29 cos p) . (33) 




FIG. 3: (Color online). Density distribution and spin texture 
of Phase II for a = —1.6. Top: Three rows are densities in xy, 
yz and xz planes respectively, three columns are for up, down 
components and the total density as explicitly labeled above 
each column. Though the total density is isotropic again, the 
density distribution for the two components exhibits more 
complex symmetry as described in the text. Bottom: 3D 
skyrmion spin texture s(r) = S(r)/n(r) of Phase II, which is 
roughly a 7r/2 rotation about y -axis of that in Phase I. The 
topological structure of the spin texture is protected by the 
time reversal symmetry. 


The average spin polarization along z axis is zero, i.e. 
(S z ) = 0. The spin texture S(r)/n(r) is presented in 
Fig. [3] The spin density in this case forms a torus 
near the yz plane and the fountain-like streamlines of 
spin pass through the hole of the torus, which is more or 
less like a 7 t/ 2 rotation of the torus in Phase I about y 
axis. Similarly, this ground state can be obtained from 
a local spin rotation from the spinor order parameter 
Cx = (c+,c_) r = ^(1,1) T that describes a system with 
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all spins pointing to the positive x direction, i.e. 

’F x = exp (-id (r) -s) ^n(r)Q x . (34) 

The spinor wavefunction Q x is related to by a 7r/2 ro¬ 
tation around y. Owing to the non-Abelian nature of 
SO(3) rotation, the spin texture of Phase II is different 
from the n/2 rotation around y of Phase I. The differ¬ 
ence between these two textures lies in the fact that the 
spin in the torus of Phase I revolves the z axis follow¬ 
ing an elliptical (oval) orbits that rotate gradually like 
the perihelion precession in celestial mechanics, while in 
Phase II the orbits are closed loops. Apart from this, 
they indeed share the same topology determined by the 
same n as can be seen from the spin streamline plot in 
Fig. El because the topological spin texture is protected 
by time reversal symmetry of the system. This skyrmion 
spin texture is proposed as the ground state in the regime 
of c > 1 in Ref. 

We find in this study that in both phase I and phase 
II, the densities of the two components are spatially sep¬ 
arated in three dimensions. We thus come up with a 
conclusion that phase separation of the spin components 
generally exists in lDj30]. 2D@|1C1 and 3D SO coupled 
Boson gases. In our case it is the SOC-induced p wave 
spatial mode involving in the variational order parameter 
that drives the two spin components spatially separated. 
As pointed out by Battye et al. [§l| phase separation 
is a prerequisite for existence of stable skyrmion, which 
explains why skyrmion spin texture appears in the varia¬ 
tional ground state of our model. Furthermore the topol¬ 
ogy of the skyrmion texture is protected by the time re¬ 
versal symmetry of the system even the phase transition 
drastically changes the density structure. 


V. SUMMARY 


We have investigated variationally the ground state 
phase diagram of weakly 3D spin-orbital coupled two- 
component Bose gases in a harmonic trap. Two phases 
for the ground state are identified depending on in¬ 
traspecies and interspecies interaction strength and the 
corresponding density distribution and the spin tex¬ 
ture are illustrated for optimized variational parameters. 
Phase I is featured with the parity symmetric and ro¬ 
tational symmetric density distribution of both spin-up 
and spin-down components and skyrmion spin texture 
with torus in the xy plane and spin streamline passing 
through the central region, while Phase II is character¬ 
ized with density distribution possessing discrete ir ro¬ 
tational symmetry about y axis and 7r rotational-time- 
reversal symmetry about x and z axis and the similar 
spin torus is in the yz plane, roughly a 7 t/ 2 rotation of 
that in Phase I about y axis. In both phases, the density 
of two components is essentially phase separated. With 
increasing interaction strength, interesting phase transi¬ 
tion occurs between the two phases, while the topology 
of ground state spin textures is protected. 
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